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SUMMARY 

Monte Carlo calculations of the self -shielding factors for multiple -body 
concentric cylindrical geometries are made and the results are compared vith 
the rational approximation and other analytic results. The blackness for a 
monoenergetic neutron flux incident upon an infinite length cylindrical shell 
of absorber vith central internal rods of absorber and multiple concentric ring 
configurations are calculated by means of the Monte Carlo method. A cosine 
incident-angular-neutron distribution and one collision model are chosen. The 
probability of absorption is calculated for thin shell and central rod combi- 
nations and for multiple concentric cylinders for a range of macroscopic ab- 
sorption cross sections. 

The self-shielding factors computed for these configurations are compared 
vith the results of the rational approximation as veil as other analytic re- 
sults . For large variations in internal structure the self -shielding factors 
shov large deviations from the rational approximation and from exact calcula- 
tions for equivalent solid rods. For multiple concentric annuli, hovever, the 
correlation of self- shielding factor vith volume -to-surface ratio is quite pre- 
cise and agrees vith the corresponding equivalent solid rod. 

To further illustrate the technique, a Monte Carlo calculation is made for 
the blackness of arrays of concentric cylinders for a single Breit Wigner reso- 
nance absorber exposed to a neutron sloving-dovn distribution. 


INTRODUCTION 

In the design of heterogeneous nuclear reactors it is essential to have 
convenient and accurate means for estimating self-shielding effects. If the 
geometry of the region for vhich the self- shielding factor is sought is moder- 
ately complex, the calculation of this parameter becomes difficult if not im- 
possible to compute by analytic means. For lumped absorbers of arbitrary shape 
and mass distribution, various approximations (refs. 1 and 2) have been used 
vith varying degrees of success . Perhaps the most videly used approximation 
has been the rational approximation of Wigner, et al. (ref. 2). The rational 
approximation asserts that, for a spatially uniform isotropic source of neutrons 



in a volume 


V, the escape probability P Q may be estimated by the formula 
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where is the total macroscopic cross section and V/S is the volume -to- 

surface ratio of the lump with 4V/S corresponding to the average chord 

length. The rational approximation, however, is subject to as much as 18 per- 
cent error (ref. l) for cylindrical rods. The error associated with the use of 
the rational approximation to determine the self -shielding of cylindrical 
shells is likewise rather large. The precision of available methods for the 
calculation of self- shielding factors for a configuration of concentric cylin- 
drical shells is unknown. The present calculations have been made by Monte 

Carlo methods to ascertain the validity of the correlation of the v/S ratio 

to the self- shielding factor for large variations in the internal structure of 
infinitely long multiple-body concentric cylindrical configurations. To fur- 
ther illustrate the utility of the Monte Carlo technique, a blackness factor 
averaged over a single-level Breit-Wigner resonance is computed for an array of 
concentric hollow cylinders. 

A Monte Carlo Fortran IV code that conveniently handles complex cylindri- 
cal geometries has been written and used in the present analysis. A non- 
scattering model has been used with a cosine angular-neutron-flux distribution 
incident upon the cylindrical arrays. 


SYMBOLS 


E energy, ev 

E . ,E minimum and maximum energy considered 

E neu energy of neutron 

F s self-shielding factor 

f optical length through cylinder 

m number of annuli 

E r r^h random number whose range is 0.0 to 1.0 

n number of neutron histories 

P probability of nontransmission 

P Q escape probability 

R distance from point to center of cylindrical assembly 
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Rj_ outside radius of shell 

Rg inside radius of shell 

R3 radius of internal rod 

R4 outside radius of single annulus (fig* 11 ) 

R5 inside radius of single annulus (fig. 11 ) 

rgi outside radius of i^* 1 annulus 

r 2i-l inside radius of i^* 1 annulus 

r 2m outside radius of annular configuration of m annuli 
S outside surface area exposed per unit length 

S a outside surface area exposed per unit length for annular configuration 

T probability of transmission 

t thickness of thin shell in thin shell and internal rod configuration 

V volume of absorbing material per unit length 

volume of absorbing material per unit length in multiple annuli configu- 
ration 

5 error 

0 angle between diameter through source and projection of neutron trajec- 

tory on plane orthogonal to axis of symmetry through source 

0 angle between diameter and tangent to inside of hollow shell (see 

1 fig. 2) 

6 U unit vector in 0 direction 

A total path length to collision 

A projection of path length to collision on cutting plane 

ir 

p mean 

s (f) z a 

2 a macroscopic absorption cross section 

Z_k total macroscopic cross section 
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standard deviation 


°s 

0 angle between neutron trajectory and projection of neutron trajectory on 

plane orthogonal to axis of symmetry through source 

unit vector in $ direction 

cp a neutron flux at surface of cylinder 

cp average flux within body 


ANALYSIS 


The Monte Carlo program assumes a nonseattering model, that is, each col- 
lision constitutes an absorption. The cross sections used are either mono- 
energetic or correspond to single-level resonances. The ratio of the number of 
absorbed neutrons to the total number of neutrons incident, which is the black- 
ness of the configuration, is calculated and appears as output from the Monte 
Carlo program; from the calculated blackness of a configuration one may calcu- 
late the self- shielding factor. 


An outline (fig. l) chronologically associated with the flow diagram of 
the Monte Carlo techniques used in this program will be followed in the analy- 
sis . A neutron is generated with 
energy taken randomly from the 
neutron energy distribution l/E 
in the interval from the minimum 
energy to the maximum energy con- 
sidered. The following rejection 
technique is used to produce this 
distribution. First a neutron 
energy is chosen in the energy in- 
terval from E m in to by the 

following formula: 



•^neu ~ ®min + ^r(®max ~ ^min) 

Next the value of N r +]_ is com- 
pared to the value of the follow- 
ing expression: 
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Figure 1. - Generalized flow chart. If the value of 

than or equal to 

expression, the value is accepted; otherwise the value E neu 


N r+ i is less 
the value of the 
is rejected and 


4 








II 


the entire process is repeated. The average of the process checked against the 
expectation value is found to give excellent agreement, better than 0.1 percent 
for ten thousand histories. 

A cosine angular distribution of the incident -neutron number flux is used 
as -well as an isotropic incident-neutron distribution. The method used bo 
generate these distributions is similar to that of Cashwell (ref. 3). The ex- 
pressions used are derived in appendix A. 

An option of the program allows a constant cross section or resonance 
parameter to be read in. When the single-level resonance absorption cross 
sections are used, Doppler broadening is neglected, and the cross sections are 
computed from the standard Breit-Wigner formula (refs. 4 and 5). In computing 
the blackness over a representative Breit-Wigner resonance, a maximum energy 
and a minimum energy are chosen. The program computes a blackness averaged 
over a l/E neutron energy spectrum. The path length to collision is calcu- 
lated in the usual manner 


A = 


In N r 
2 t 


The thickness of material seen by the neutron (the optical thickness) f is 
then calculated. This thickness of material is compared with the length A. 

If A is longer than f, the neutron is not absorbed. If A is less than or 
equal to f, the neutron is absorbed. The value of f for the hollow cylinder 
is computed from 

2 

f = - 

cos 0 

and for the thin shell and rod (fig. 2) from 

2^ cos 9 - (R| - Rj sin 2 0) ly/2 + (R§ - r| sin 2 ©) 1 / 2 ] 

COS $ ~ ~ 


0 - (Bg - Bp sin 0) 


The derivation of these expressions is shown in appendix B. 

A self -shielding factor F s is defined to be the ratio of the average 
flux inside a body to the surface flux (ref. 6) 


The inwardly directed current or number of neutrons per second striking a body 
immersed in an isotropic scattering medium is cp s S/4, where S is the outside 
surface area of the body containing material of volume V and of macroscopic 
absorption cross section £ a . If T is the fraction of the incident neutrons 
leaving the assembly, the blackness 1 - T is given by 
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or 


where 




1 - T _ 1 - T 

, 4V* 2| 
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‘ Configurations 

Hollow cylinder and shell and rod configuration . - Self-shielding calcu- 
lations are made for two geometric configurations. The thin cylinder with an 
internal rod shown in figure 2 and the hollow cylinder are chosen because they 

represent two extremes of internal structure. 
Since these configurations are not completely 
specified by the parameter a particular 
family is chosen to illustrate the differences 
in self- shielding factors . The surface area 
exposed per unit length is the same for both 
configurations. The internal volume and den- 
sity of absorber per unit length of both as- 
semblies is held constant for each value of | . 

For a hollow cylinder of outside radius 
Rp, it is shown in reference 7 that 

i = ^l^a cos2 ^l 

where 6^ is angle between a diameter and 
a line tangent to the inside wall of the hol- 
low cylinder (see fig. 2). Thus, to vary 
one can vary either the outside radius , the 
macroscopic absorption cross section, or the 
amount of material inside the shell. The hol- 
low cylinders calculated used the product 
Figure 2. - Thin cylinder with internal rod configuration. = 1*0, a constant; thus, to vary |, the 

amount of material inside the hollow shell is 
changed, that is, the central hollow gradually is filled with material. As 
0 1 o, | -* 1.0} for £ = 1.0, under the foregoing conditions, a solid rod is 
obtained. 

The thin shell and rod cases shown in figure 2 used constant values of the 
ratio of shell thickness to outside radius t/Rp of 0.0005 and 0.05. The dif- 
ference in mass between the thin shell and the hollow cylinder configuration 
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appeared as the central rod of the thin shell and the central rod configura- 
tion. Representative values of dimensions are shown in table I. 


Multiple-body concentric cylindrical configurations . - Self-shielding cal- 
culations are made for a series of concentric shells illustrated in figure 3. 

The number of shells making 
TABLE I. - REPRESENTATIVE DIMENSIONS FOR HOLLOW CYLINDERS up thg assem bly is a param- 

and thin shells with internal rods eter varied from 1 to 20; 

the ratio of shell material 


Hollow cylinder 

Thin shell and 





internal rod 

Surface- 

Outside 

Radius of 

Ratio of shell 

Radius 

to -volume 

radius 

inner hoi- 

thickness to 

of in- 

ratio of 

of hollow 

low shell,* 

outside 

radius , 

ternal 

lump. 

cylinder. 

r 2> 

t/Rn 

rod, 

s/v. 

R i> 

cm 




cm”-*- 

cm 




cm 

2000 

1.000 

0.9995 

0.0005 


1667 



.9994 



0.01414 

1000 



.9990 



.03161 

500 



.9980 



.05474 

250 



.9960 



.08357 

100 



.9900 



.13748 

50 



.9800 



.19647 

25 



.9600 



.27821 


volume to outside surface 
v/s is held constant. The 
volume per unit length of 
the annuli is given by 

m 

V A = n ^2 ( r 2i + r 2i-l) 
i=l 


x ( r 2i " r 2i-l) 

where m is the number of 
annuli, r 2i is the outside 
radius of the i^* 1 shell, 
and rg i_i is the inside 
radius of the i^* 1 shell. 


The outside surface area per unit length is 



For a given | there is a rather varied choice of geometries; that is, | does 
not completely specify an annular cylindrical geometry. The following ancil- 
lary conditions were added to specify the geometry of the annuli. 
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(1) The cylindrical thickness for a given configuration containing m 
cylinders shall be equal; that is, 

rgq - r 2i-l = "k (i = 1,2,3,. . .,m) 

(2) The spacing between cylinders for a given m shall be equal. For 
example, the distance from the z-axis of symmetry to the outside of ring 1 is 
just one-half that to the outside of ring 2. 

(3) The ratio V/S is held constant, while the number of rings is allowed 
to vary. 

( 4 ) For a given m and V/S, 2 a is varied to change the value of g. 


( 5 ) The outside radius of all concentric cylindrical configurations con- 
sidered herein is chosen at 1.5 inches. 


TABLE II. - DIMENSIONS OB MULTIPLE 


CONCENTRIC CYLINDERS 


Number 

Outside ra- 

Thickness 

Spacing of 

of 

dius of con- 

of annuli. 

annuli (out 

rings 

figuration , 

cm 

side radius 


cm 


to outside 





radius ) , 





cm 

20 

3. 

81 

0.0240 

0.1905 

10 



.0461 

.3810 

6 



.0728 

.6350 

4 



.1024 

.9525 

3 



.1280 

1.2700 

2 



.1720 

1.9050 

1 



.2590 



are shown in table IF. To illustrate 
inc ident -angular -flux distributions , 


The foregoing conditions spec- 
ify the concentric ring geometries. 
The values of the radii chosen are 
shown in table II. The values of 
m, g, and F s are shown in 
table III. A cosine incident angu- 
lar neutron number flux is used. 


Concentric Cylinders of 

Resonance Absorbers 

The Monte Carlo program is also 
executed for a representative reso- 
nance (AglOS at 5.20 ev) (see 
ref. 4). Values of the radii chosen 
the sensitivity of the results to 
an isotropic incident-neutron case was 
also calculated. 


TABLE III. - SELF-SHIELDING FACTORS FOR 


Statistical Analysis 


CYLINDRICAL CONFIGURATIONS 


Number 


Values of 

6 a 


of 

rings 

0.10 

0.20 

0.30 

0.40 

0.50 

1 

0.875(0) 

0.775(2) 

0.696(3) 

0.631(6) 

0.574(0) 

2 

.885(5) 

.789(2) 

.708(6) 

.642(0) 

.584(9) 

3 

.900(7) 

.796(9) 

.713(3) 

.645(7) 

.587(4) 

4 

.891(5) 

.796(1) 

.720(4) 

.652(1) 

.589(8) 

6 

.891(4) 

.797(2) 

.719(1) 

.651(1) 

.594(5) 

10 

.894(5) 

.804(2) 

.721(8) 

.652(0) 

.593(9) 

20 

.902(5) 

.802(9) 

.720(3) 

.650(0) 

.594(2) 

00 

.88502 

.79303 

. 71649 

.65162 

.59595 


a Cosine incident. 


Given a neutron, there is a 
probability T that it shall be 
transmitted through the assembly 
and a probability P = 1 - T 
that it will not be transmitted. 
Let n stochastically independ- 
ent neutrons be injected into 
the system; then, x neutrons 
will be transmitted. These are 
precisely the postulates neces- 
sary for the generation of a 
binomial distribution (ref. 8). 


8 




TABLE IV. - DIMENSIONS OF -RESONANCE CONFIGURATION 


Number 

Outside ra- 

Thickness 

Spacing of 

of 

dius of con- 

of annuli , 

annuli ( out- 

rings 

figuration, 

cm 

side radius 


cm 


to outside 





radius ) , 





cm 

20 

1.000 

0.0005236 

0.0500000 

18 



.0005787 

. 0555556 

16 



.0006469 

.0625000 

14 



. 0007332 

. 0714286 

12 



.0008461 

. 0833333 

10 



.0010000 * 

. 1000000 

8 



.0012224 

.1250000 

6 



.0015721 

.1666667 

5 



.0018345 

.2000000 

4 



.0022019 

.2500000 

3 



.0027532 

.3333333 

2 



.0036723 

.5000000 

1 



. 0055102 



In order to calculate the first 
two moments of the binomial distribu- 
tion, the following formulas are used: 

[x = nT 

where \x is the mean and 
a s = -|/nTP 

where a s is the standard deviation. 
Thus, the best estimate of T is 

T = — 
n 

and the best estimate of o s is 
o s = "v^lP 


Hence, the fractional error associated with one standard deviation is 


Since 



a s 

The fractional error of F s is 

5F s 1 5T 

F s 4V T 

s 

For example, if n = 10,000, T = 0.75, and [s a ( 4V r /s)] _ ‘ L = 3.33, 

5F s _ / 47\ _1 /T 

F s " V a SJ fnT 


= 3.33 X 0.57 percent 


= 1.90 percent 


Hence 1.90 percent is the percentage error associated with one stan dar d devia- 
tion. The errors given to the Monte Carlo results have been calculated in this 
manner and are standard deviations . These standard deviations are somewhat 
smaller than the symbols marking the data unless otherwise shown. 
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RESULTS AND DISCUSSION 


Self-shielding factors calculated for two extreme geometric configura- 
tions , the thin cylinder with an internal rod and the hollow cylinder (see 
fig. 2, p. 6) are presented in figure 4 as a function of A large differ- 
ence in F s for these cases is obtained at the smaller values of These 



calculated values are compared with exact results for rods (refs. 7 and 9). 

For the chosen geometries, the configurations correspond to rods at | = 1 and 
cross the exact solution here. 

The hollow cylinder has its mass concentrated toward its periphery, while 
the thin cylinder and internal rod has its mass concentrated near the center of 
the structure; for a cosine incident -angular -flux distribution, therefore, the 
average chord length is longer for the latter case, and thus the self- shielding 
effect is greater. (As can be seen from figure 4, the largest discrepancy in 
the correlation lies in the range 0.00 < | ^ 0.50.) It is evident that the 
parameter £ = (2V/s)Z a does not correlate the self -shielding factors in cyl- 
indrical configurations with large variations of internal structure. 

For multiple body concentric cylinders, results are shown in figures 5 
and 6 for the range £ = 0.00 to £ = 0.50. As noted before, the annuli are 
chosen all of the same thickness and equally spaced. The outside diameter, a 
constant, chosen for the outermost annulus was 3.0 inches. The annular dimen- 
sion remained a constant (see table II, p. 8), and the cross section £ a was 
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— Value of self-shieidmg for 
rods of equal ^ 


Number of annuli, m 

Figure 5. - Self-shielding factors plotted against number of annuli. 

varied to change g. Results for F s are shown in table III. As can be seen 
in figure 5, the self -shielding factors tend to approach the solid rod asymp- 
tote for large values of £. For low values of £, the self- shielding curves 
for multiple annuli exceed the solid rod asymptote. This effect may well be 
due to streaming across the internal voids. As the blackness of the array is 
increased, this effect ceases to be important. For large numbers of rings, the 
curves approach the self-shielding factors for rods . It appears that, for the 
range of variables studied, the self -shielding factors for the multiple concen- 
tric cylinders is reasonably well approximated by the equivalent solid rod 
self- shielding factors as prescribed by the £ correlation. 


The results shown in figure 7 present the blackness factors averaged over 
a Breit-Wigner resonance and a l/E neutron energy distribution. The value of 
V/S is kept constant at 0.0055 centimeter and the number of cylinders is 
varied from 1 to 20. 



Blackness factor Self-shielding factor, 



Figure 6. - Self-shielding factors for multiple cylindrical configurations. 


Blackness factors, 
rather than self-shielding 
factors, are plotted inasmuch 
as a simple representative 
value of | over an entire 
resonance cannot be obtained. 

The blackness for an 
isotropic angular-flux distri- 
bution is compared with the 
blackness obtained using a 
cosine incident angular dis- 
tribution. The isotropic 
distribution is introduced as 
a mathematical device to as- 
certain the effect of a change 
in angular distribution on 
the blackness factor averaged 
over a Breit-Wigner resonance 
line shape and a l/E energy 
spectrum. 

The change in blackness 
is more sensitive to an iso- 
tropic incident flux than to 
a cosine incident flux, as 
shown in figure 7, for smaller 
numbers of annuli. The 
greatest difference between 
the cosine angular and the 
isotropic angular -incident- 
flux distributions occurs for 
a single annulus, but the so- 




Figure 7. - Blackness factor plotted against number of annuli 
over a resonance. 


Figure 8. - Blackness factors plotted against number of annuli 
for monoenergetic flux. 
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lution appears to converge for large numbers of concentric cylinders. In order 
to ascertain whether this apparent convergence is real, monoenergetic blackness 
factors are computed and are shown in figure 8. These results indicate that, 
if the resonance calculations had been done for a larger number of annuli, the 
blackness factor curves would cross. 


Lewis Research Center 

Rational Aeronautics and Space Administration 
Cleveland, Ohio, July 29, 1964 
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APPENDIX A 


DERIVATION OF ANGULAR DISTRIBUTIONS 

The results of the following derivation were used to obtain an isotropic 
neutron-angular-flux distribution. A differential element of area dA in 
spherical coordinates can be written 

dA = r^ sin $ d4> d8 

If r = 1 and polar coordinates are used, the following identities may be in- 
troduced (see fig. 9): 


u = cos a 
v = cos p 
w = cos T 




(a) Polar coordinate system. (b) Cartesian coordinate system. 

Figure 9. - Coordinate systems. 


By simple transformation, 
r = P = 1 d$ = dy 

r = 4) d0 = dv 

v = cos p = 6 dw =* -sin 4> d4> 

w = cos $ 
and 


dA = sin 4> d4> d0 = -dw dv 

The probability density function of w, p(w) dw can be written 

r i n 2itr 2 sin 4> d4> 1 . * 1 , 

p(w; dw = ^ — sin $ d4> = — dw 


(Al) 


If N r is a random number, the solution of equation (Al) gives the method 
of choosing w for an isotropic distribution 


N r = 


i 

• / -i 


i . . i , 

- dw = 2 W 


w 


-1 


2 + V 


and 


W S= 2Ny» - 1 
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Thus w of the set u, v, and w can he selected. 

Since "$ u is orthogonal to the vector 0 U , 0 can be obtained indepen- 
dently. The fact that 6 may range from -it to jt with equiprobability 
means the value of 6 can be determined in the following manner: 

Hr - r P(e')M' - f° fr - sr 

•'-Jt *'-JT 


and 


e = rt(2W r - l) 

A cosine distribution may be gained by the following substitution: 

p(w) = 2w 

by definition and 

p(w) dw = 2w dw 

Let 




= / 


2w dw = w^ 


Hence 



APPENDIX B 


DERIVATION OF NEUTRON PATH IFNG-TH IN MATERIAL 


The equation used to calculate the path length in material shall he de- 
rived in this appendix. First consider the section at a cutting plane perpen- 
dicular to the infinite hollow rod (see fig. 10). This figure is composed of 

two concentric circles A and B with radii Rq and 
R2, respectively, as well as the projection into the 
cutting plane of a neutron trajectory A'A". 

The distance A T A n must first be obtained, and 
then the distance B'B" subtracted. Choosing the point 
A 1 as the origin of the coordinate system and writing 
the equation for A T A ,r in terms of Rq and Q result 
in the following equation: 

A*A M = 2Rq cos G (Bl) 

The equation for the circle B in terms of Rq, R^, 

and G is 



Rg = r 2 + R^ - 2rR^ cos G 


(B2) 


where (r,0) is the polar vector that describes the circle B f from A 1 , the 
origin. Solving for r gives 


r = R^ cos 


e ± y^ 2 cos 2 e - (r| - Rp 


Since cos 2 © - 1 = -sin 2 ©. 


r = R]_ cos © ± y] 

which describes the circle B. Hence 


2 2 2 
R2 - Rq sin © 


(B3) 


(B4) 


B 


•B" = 2-^i 


'R 2 - R 2 s in 2 0 


(B5) 


Thus the expression for the projection of the path length on the cutting plane 
is 


2R-^ cos G 


2 



- R^ sin 


(B6) 


When the projection intersects an inner rod or annulus, the following expres- 
sions (which may be obtained from eq . (B5)) are used (see fig. ll): 


C’C” = 2 



sin 2 © 


(B7) 
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D'D" - E’E" = 2 yR| - R^ sin^Q - yR~ - R^ sin^e \,B8 ) 

Once the length of the entire projection Ap (on the material in the cut- 
ting plane) has been computed, the chord length A is calculated by the follow 
ing formula: 

A 

cos $ 

where $ is the angle the neutron velocity vector makes with the cutting plane 



Figure 11. - Projection of neutron trajectories through various configurations. 
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